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ABSTRACT 



We discuss the late time evolution of the gravitational clustering in an 
expanding universe, based on the nonlinear scaling relations (NSR) which 
connect the nonlinear and linear two point correlation functions. The existence 
of critical indices for the NSR suggests that the evolution may proceed towards 
a universal profile which does not change its shape at late times. We begin 
by clarifying the relation between the density profiles of the individual halos 
and the slope of the correlation function and discuss the conditions under 
which the slopes of the correlation function at the extreme nonlinear end 
can be independent of the initial power spectrum. If the evolution should 
lead to a profile which preserves the shape at late times, then the correlation 
function should grow as [in a Q = 1 universe] even at nonhnear scales. We 
prove that such exact solutions do not exist; however, there exists a class of 
solutions ( "psuedo-hnear profiles", PLP's for short) which evolve as to a good 
approximation. It turns out that the PLP's are the correlation functions which 
arise if the individual halos are assumed to be isothermal spheres. They are also 
configurations of mass in which the nonlinear effects of gravitational clustering 
is a minimum and hence can act as building blocks of the nonlinear universe. 
We discuss the imphcations of this result. 

Subject headings: cosmology: theory - structure formation - two 

point correlation functions - power spectra 
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1. Introduction 

The evolution of large number of particles under their mutual gravitational 
influence is a well-defined mathematical problem. If such a system occupies 
a finite region of phase space at an initial instant, and evolves via newtonian 
gravity, then it does not reach any sensible 'equilibrium' state. The core region 
of the system will keep on shrinking and will be eventually be dominated by a 
few hard binaries. Rest of the particles will evaporate away to large distances, 
gaining kinetic energy from the shrinking core [for a discussion of such systems, 
Padmanabhan 1990| . 



see 



The situation is drastically different in the presence of an expanding 
background universe characterised by an expansion factor a{t). Firstly, the 
expansion tends to keep particles apart thereby exerting a civilising influence 
against newtonian attraction. Secondly, it is now possible to consider an 
infinite region of space filled with particles. The average density of particles 
will contribute to the expansion of the background universe and the deviations 
from the uniformity will lead to clustering. Particles evaporating from a 
local overdense cluster cannot escape to "large distances" but necessarily will 
encounter other deep potential wells. Naively, one would expect the local 
overdense regions to eventually form gravitationally bound objects, with a 
hotter distribution of particles hovering uniformly all over. As the background 
expands, the velocity dispersion of the second component will keep decreasing 
and they will be captured by the deeper potential wells. Meanwhile, the 
clustered component will also evolve dynamically and participate in, e.g mergers. 
If the background expansion and the initial conditions have no length scale, 
then it is likely that the clustering will continue in a hierarchical manner ad 
infinitum. 
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Most of the practising cosmologists will broadly agree with the above picture 
of gravitational clustering in an expanding universe. It is, however, not easy to 
translate these concepts into a well-defined mathematical formalism and provide 
a more quantitative description of the gravitational clustering. One of the key 
questions regarding this system which needs to be addressed is the following: 
Can one make any general statements about the very late stage evolution of the 
clustering ? For example, does the power spectrum at late times 'remember' the 
initial power spectrum or does it possess some universal characteristics which 
are reasonably independent of initial conditions ? [This question is closely 
related to the issue of whether gravitational clustering leads to density profiles 
which are universal. [Navarro, Frenk fc White 1996| . 



We address some aspects of this issue in this paper and show that it is 
possible to provide (at least partial) answers to these questions based on a 
simple paradigm. The key assumption we shall make is the following: Let ratio 
between mean relative pair velocity v{a,x) and the negative hubble velocity 
{—ax) be denoted by h{a,x) and let ^{a,x) be the mean correlation function 
averaged over a sphere of radius x. We shall assume that h{a, x) depends 
on a and x only through ^{a,x); that is, h{a,x) = h[C,{a,x)]. With such a 
minimal assumption, we will be able to obtain several conclusions regarding the 
evolution of power spectrum in the universe. Such an assumption was originally 
introduced — in a different form — by Hamilton ([Hamilton et al. 1991|) . The 
present form, as well as its theoretical implications were discussed in [Nityananda] 



fc Padmanabhan 1994[ , and a theoretical model for the scaling was attempted by 
Padmanabhan ( [Padmanabhan 1996[ ). It must be noted that simulations indicate 
a dependence of the relation h{a,x) = h[C,{a,x)] on the intial spectrum and also 
on cosmological parameters [[Peacock fc Dodds 1994[ ; [Peacock fc Dodds 1996[ ; 
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Padmanabhan et al. 199"6| ; |Mo, Jain fc White 1995 |. Most of our discussion 



in this paper is independent of this fact or can be easily generahsed to such 
cases. When we need to use an exphcit form for h we shall use the original ones 
suggested by Hamilton ( [Hamilton et al. 1991|) because of its simplicity. 

Since this paper addresses several independent but related questions, we 
provide here a brief summary of how it is organised. Section 2 studies some 
aspects of nonlinear evolution based on the assumption mentioned above. We 
begin by summarising some previously known results in subsection 2.1 to set 
up the notation and collect together in one place the formulas we need later. 
Subsection 2.2 makes a brief comment about the critical indices in gravitational 
dynamics so as to motivate later discussion. In section 3, we discuss the relation 
between density profiles of halos and correlation functions and derive the 
conditions under which one may expect universal density profiles in gravitational 
clustering. In section 4 we show that gravitational clustering does not admit self 
similar evolution except in a very special case. We also discuss the conditions for 
approximate self-similarity to hold. Section 5 discusses the question as whether 
one can expect to find power spectra which evolve preserving their shape, 
even in the nonlinear regime. We first show, based on the results of section 4, 
that such exact solutions cannot exist. We then discuss the conditions for the 
existence of some approximate solutions. We obtain one prototype approximate 
solution and discuss its properties. The solution also allows us to understand 
the connection between statistical mechanics of gravitating systems in the small 
scale and evolution of correlation functions on the large scale. Finally, section 6 
discusses the results. 
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2. General features of nonlinear evolution 

Consider the evolution of the system starting from a gaussian initial 
fluctuations with an initial power spectrum, Pin{k). The fourier transform of 
the power spectrum defines the correlation function ^ (a, x) where a oc t^/^ is the 
expansion factor in a universe with = 1. It is more convenient to work with 
the average correlation function inside a sphere of radius x, defined by 

3 /"^ 

^ia,x) = ^ i{a,y)y'^dy (1) 
This quantity is related to the power spectrum P(a, k) by 

^(x, a) = / —P{a,k)[sm{kx) — kxcos{kx)] (2) 

ZTT^X'^ JO k 

with the inverse relation 

An f°° 

P(a,k) = — / dx X ^{a, x) [sm{k x) — k X cos{k x)] (3) 
3k JO 

In the linear regime we have ^l{o,,x) oc a^^j„(ai,x). 

We now recall that the conservation of pairs of particles gives an exact 
equation satisfied by the correlation function ([Peebles 19801 ): 



where v{a,x) denotes the mean relative velocity of pairs at separation x and 
epoch a. Using the mean correlation function ^ and a dimensionless pair velocity 
h{a,x) = —{v/dx), equation (H) can be written as 

' ^ h-^) (1 + = 3^1 + (5) 



9 In a 9 In X 
This equation can be simplified by introducing the variables 

v4 = lna, X = lna;, D{X, A) = \n{l + ^) (6) 
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in terms of which we have ( Nityananda fc Padmanabhan 199i ) 



dD dD 

^-MAX)— = 3MAX) (7) 

At this stage we shall introduce our key assumption, viz. that h depends on 
{A^X) only through ^ (or, equivalently, D). Given this single assumption, 
several results follow which we shall now summarise. 



2.1. Formal solution 



Given that h = h[C,{a, x)], one can easily integrate the equation (|^) to find the 
general solution [see [Nityananda fc Padmanabhan 199^ ] . The characteristics of 



this equation (|^) satisfy the condition 

x'{l + = l' (8) 

where / is another length scale. When the evolution is linear at all the relevant 
scales, ^ ^ 1 and I x. As clustering develops, ^ increases and x becomes 
considerably smaller than /. The behaviour of clustering at some scale x is 
then determined by the original linear power spectrum at the scale / through 
the "flow of information" along the characteristics. This suggests that we can 
express the true correlation function ^{a, x) in terms of the linear correlation 
function evaluated at a different point. This is indeed true and the 

general solution can be expressed as a nonlinear scaling relation (NSR, for short) 
between ^^(a, /) and ^{a,x) with / and x related by equation(H). To express 
this solution we define two functions V{z) and U{z) where V{z) is related to the 
function h{z) by 



- 8 - 



and U{z) is the inverse function of V{z). Then the solution to the equation 
can be written in either of two equivalent forms as: 



(10) 



where 



X 



[1 + C,) ( [Nityananda fc Padmanabhan 199^ . Given the form of 



h{C,) this allows one to relate the nonlinear correlation function to the linear one. 

From general theoretical considerations [see Padmanabhan 1996|| it can be 
shown that V{z) has the form: 



Viz) 



111 



1 < 1) 

^1/3 [1 < z < 200) ( 
;z2/3 (200 < z) 

In these three regions h{z) ^ [(22;/3), 2, 1] respectively. We shall call these 
regimes, linear, intermediate and nonlinear respectively. More exact fitting 
functions to V{z) and U{z) have been suggested in literature, [see [Hamilton et 



al. 19911 ; |Mo, Jain fc White 1995] [Peacock fc Dodds 1994| ]. When needed in this 

paper, we shall use the one given in Hamilton et al.,1991: 

/ 1 + 0.0158 22 + 0.000115^3 



1/3 



V{z) 



Uiz) 



0.926 ^2 - 0.0743 z^ + 0.0156 z^ ^ 
2 + 0.358 2^ + 0.0236 z^ 



1 + 0.0134 23 + 0.0020 2^/2 
Equations (llOD and ([T^Jl3|) implicitly determine .^(a,a;) in terms of ^^(a,^) 



(12) 
(13) 



2.2. Critical indices 

These NSR already allow one to obtain some general conclusions regarding 
the evolution. To do this most effectively, let us define a local index for rate of 
clustering by 

(14) 

dm a 
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which measures how fast ^(a,a;) is growing. When -C 1, then = 2 

irrespective of the spatial variation of ^ (a, x) and the evolution preserves the 
shape of ^{a,x). However, as clustering develops, the growth rate will depend 
on the spatial variation of ^(a, a;). Defining the effective spatial slope by 

-K„K^)+31.^^i^ (15) 

one can rewrite the equation (|^) as 

= - n.ff) (16) 

At any given scale of nonlinearity, decided by C,{a,x), there exists a critical 
spatial slope tIc such that Ua > 2 for Ueff < ric [implying rate of growth is 
faster than predicted by linear theory] and < 2 for rig// > [with the rate 
of growth being slower than predicted by linear theory]. The critical index is 
fixed by setting na = 2 in equation ( ]T^ ) at any instant. This feature will tend 
to "straighten out" correlation functions towards the critical slope. [We are 
assuming that ^(a, x) has a slope that is decreasing with scale, which is true for 
any physically interesting case]. From the fitting function it is easy to see that 
in the range 1 ^ ^ 200, the critical index is nc ~ —1 and for 200 ^ C,, the 
critical index is ric ^ —2 ( |Bagla fc Padmanabhan 1997| ). This clearly suggests 



that the local effect of evolution is to drive the correlation function to have a 
shape with (l/x) behaviour at nonlinear regime and (l/x^) in the intermediate 
regime. Such a correlation function will have ~ 2 and hence will grow at a 
rate close to a^. We shall say more about this in section 3 below. 

3. Correlation functions, density profiles and stable clustering 

Now that we have a NSR giving ^(a, x) in terms of ^^(a,/) we can ask the 
question: How does f (a, x) behave at highly nonlinear scales or, equivalently, at 
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any given scale at large a ? 

To begin with, it is easy to see that we must have v = —ax or h = 1 for 
sufficiently large ^{a, x) if we assume that the evolution gets frozen in proper 
coordinates at highly nonlinear scales. Integrating equation (^) with h = 1, we 
get ^{a,x) = a^F{ax); we shall call this phenomenon "stable clustering". There 
are two points which need to be emphasised about stable clustering: 

(1) At present, there exists some evidence from simulations ( [Fadmana bhan| 
et al. 1996 ) that stable clustering does not occur in a Q = 1 model. In a formal 
sense, numerical simulations cannot disprove [or even prove, strictly speaking] 
the occurrence of stable clustering, because of the finite dynamic range of any 
simulation. 

(2) . Theoretically speaking, the "naturalness" of stable clustering is often 
overstated. The usual argument is based on the assumption that at very small 
scales — corresponding to high nonlinearities — the structures are "expected 
to be" frozen at the proper coordinates. However, this argument does not take 
into account the fact that mergers are not negligible at any scale in an = 1 
universe. In fact, stable clustering is more likely to be valid in models with i7 < 1 
— a claim which seems to be again supported by simulations ( Padmanabhan et 



al. 199q ). 



// stable clustering is valid, then the late time behaviour of ^ (a, x) cannot 
be independent of initial conditions. In other words the two requirements: (i) 
validity of stable clustering at highly nonlinear scales and (ii) the independence 
of late time behaviour from initial conditions, are mutually exclusive. This is 
most easily seen for initial power spectra which are scale- free. If Pin{k) oc A;" so 
that ^^(a, x) oc a'^x~^"'~^^\ then it is easy to show that ^(a,x) at small scales will 
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vary as 



- 6 3{n+3) - 

^{a,x) (X a'^x~~^] (C> 200) (17) 



if stable clustering is true. Clearly, the power law index in the nonlinear regime 
"remembers" the initial index. The same result holds for more general initial 
conditions. 

What does this result imply for the profiles of individual halos? To answer 
this question, let us start with the simple assumption that the density field 
p(a, x) at late stages can be expressed as a superposition of several halos, each 
with some density profile; that is, we take 

P(«>x) = ^/(x-Xi,a) (18) 

i 

where the i-th halo is centered at Xj and contributes an amount /(x — Xj,a) 
at the location Xj [We can easily generalise this equation to the situation in 
which there are halos with different properties, like core radius, mass etc by 
summing over the number density of objects with particular properties; we 
shall not bother to do this. At the other extreme, the exact description merely 
corresponds to taking the /'s to be Dirac delta functions]. The power spectrum 
for the density contrast, (5(a,x) = {p/pb — 1), corresponding to the p(a, x) in 



18]) can be expressed as 



P(k,a) cx {a'\f{k, 



a) 



(a3|/(k,a)|)'Pcent(k,a) (20) 



exp — ik ■ Xj(a) 



(19) 



oc 



where Pcent(k, a) denotes the power spectrum of the distribution of centers of 
the halos. 

If stable clustering is valid, then the density profiles of halos are frozen 
in proper coordinates and we will have /(x — Xj,a) = f{a (x — Xj)); hence 
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the fourier transform will have the form /(k, a) = /(k/a). On the other 
hand, the power spectrum at scales which participate in stable clustering must 
satisfy P(k, a) = P(k/a) [This is merely the requirement C,{a,x) = a^F{ax) 
re-expressed in fourier space]. From equation ( ^OD it follows that we must have 
-Pcent(k, a) = constant independent of k and a at small length scales. This can 
arise in the special case of random distribution of centers or — more importantly 
— because we are essentially probing the interior of a single halo at sufficiently 
small scales. [Note that we must necessarily have Pcent ~ constant, for length 
scales smaller than typical halo size, by definition]. We can relate the halo 
profile to the correlation function using pO|). In particular, if the halo profile 
is a power law with / oc r~^, it follows that the .^(a, x) scales as x~'^ [ see also 
McClelland fc Silk 1977| ; |Sheth fc Jain 1997|| where 



7 = 2e-3 (21) 

Now if the correlation function scales as [— 3(n + 3)/(n + 5)], then we see 
that the halo density profiles should be related to the initial power law index 
through the relation 

. = 5(!i±i) (22) 

So clearly, the halos of highly virialised systems still "remember" the initial 
power spectrum. 

Alternatively, one can try to "reason out" the profiles of the individual halos 
and use it to obtain the scaling relation for correlation functions. One of the 
favourite arguments used by cosmologists to obtain such a "reasonable" halo 
profile is based on spherical, scale invariant, collapse. It turns out that one 
can provide a series of arguments, based on spherical collapse, to show that — 
under certain circumstances — the density profiles at the nonlinear end scale as 
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[— 3(?2 + 3)/(n + 5)] . The simplest variant of this argument runs as follows: If we 
start with an initial density profile which is r""", then scale invariant spherical 
collapse will lead to a profile which goes as with /5 = 3a/(l + a) [see eg., 
Padmanabhan, 1996, 1996a and references cited therein]. Taking the intial slope 
as a = (ra + 3)/2 will immediately give /? = 3{n + 3)/(n + 5). [Our definition of 
the stable clustering in the last section is based on the scaling of the correlation 
function and gave the slope of [— 3(?2 + 3)/{n + 5)] for the correlation function. 
The spherical collapse gives the same slope for halo profiles.] In this case, when 
the halos have the slope of e = 3(?7. + 3)/{n + 5), then the correlation function 
should have slope 

(23) 

n + 5 

Once again, the final state "remembers" the initial index n. 

Is this conclusion true ? Unfortunately, simulations do not have sufficient 
dynamic range to provide a clear answer but there are some claims [see [Navarro,] 
Frenk fc White 1996 ] that the halo profiles are "universal" and independent 
of initial conditions. The theoretical arguments given above are also far from 
rigourous (in spite of the popularity they seem to enjoy!). The argument for 
correlation function to scale as [— 3(?2 + 3)/(n + 5)] is based on the assumption of 
h = 1 asymptotically, which may not be true. The argument, leading to density 
profiles scaling as [— 3(r;, + 3)/{n + 5)], is based on scale invariant spherical 
collapse which does not do justice to nonradial motions. Just to illustrate the 
situations in which one may obtain final configurations which are independent 
of initial index n, we shall discuss two possibilities: 

(i) As a first example we will try to see when the slope of the correlation 
function is universal and obtain the slope of halos in the nonlinear limit using 



our relation (|2T|) . Such an interesting situation can develop if we assume that 
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h reaches a constant value asymptotically which is not necessarily unity. In 
that case, we can integrate our equation (|^) to get ^(a, x) = a^^F[a^x] where 
h now denotes the constant asymptotic value of of the function. For an initial 
spectrum which is scale-free power law with index n, this result translates to 



where 7 is given by 



— 27 

C{a,x) oc a^x""^ (24) 



_ 3h{n + 3) 
^-2 + h{n + 3) ^^^^ 



We now notice that one can obtain a 7 which is independent of initial power law 
index provided h satisfies the condition h{n + 3) = c, a constant. In this case, 
the nonlinear correlation function will be given by 



3c 



^(a,a;) OC a (2+^^T(^a; 2+'= (26) 
The halo index will be independent of n and will be given by 

Note that we are now demanding the asymptotic value of h to explicitly depend 
on the initial conditions though the spatial dependence of ^(a,x) does not. 
In other words, the velocity distribution — which is related to h — still 
"remembers" the initial conditions. This is indirectly reflected in the fact that 
the growth of ^{a,x) — represented by a^^^^^'^^'^^^'^'^^'^^ — does depend on the 
index n. 

As an example of the power of such a — seemingly simple — analysis note 
the following: Since c > 0, it follows that e > (3/2); invariant profiles with 
shallower indices (for e.g with e = 1) are not consistent with the evolution 
described above. 



- 15 - 



(ii) For our second example, we shall make an ansatz for the halo profile and 
use it to determine the correlation function. We assume, based on small scale 
dynamics, that the density profiles of individual halos should resemble that of 
isothermal spheres, with e = 2, irrespective of initial conditions. Converting this 
halo profile to correlation function in the nonlinear regime is straightforward 
and is based on equation (PT]): If e = 2, we must have 7 = 2e — 3 = 1 at small 
scales; that is ^(a, x) oc at the nonlinear regime. Note that this corresponds 
to the critical index at the nonlinear end, ne// = ric = —2 for which the growth 
rate is — same as in linear theory. [This is, however, possible for initial power 
law spectra, only if e = 1, i.e h{n + 3) = 1 at very nonlinear scales. Testing the 
conjecture that h{n + 3) is a constant is probably a little easier than looking for 
invariant profiles in the simulations but the results are still uncertain]. 

The corresponding analysis for the intermediate regime, with 
1 ^ ^{a,x) ^ 200, is more involved. This is clearly seen in equation 
( POD which shows that the power spectrum [and hence the correlation function] 
depends both on the fourier transform of the halo profiles as well as the power 
spectrum of the distribution of halo centres. In general, both quantities will 
evolve with time and we cannot ignore the effect of Pcent(^, cl) and relate P{k, a) 
to f{k,a). The density profile around a local maxima will scale approximately 
as p oc ^ while the density profile around a randomly chosen point will scale as 
p oc [The relation 7 = 2e — 3 expresses the latter scaling of ,^ oc p^]. There 
is, however, reason to believe that the intermediate regime (with 1 ^ ^ 200) 
is dominated by the collapse of high peaks ( [Padmanabhan 1996[ ) . In that case, 
we expect the correlation function and the density profile to have the same slope 
in the intermediate regime with ^(a, x) oc (1/x^). Remarkably enough, this 
corresponds to the critical index ng// = = —1 for the intermediate regime for 
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which the growth is proportional to a? . 

We thus see that if: (i) the individual halos are isothermal spheres with 
(1/x^) profile and (ii) if oc p in the intermediate regime and ^ oc in the 
nonlinear regime, we end up with a correlation function which grows as a? at all 
scales. Such an evolution, of course, preserves the shape and is a good candidate 
for the late stage evolution of the clustering. 

While the above arguments are suggestive, they are far from conclusive. It 
is, however, clear from the above analysis and it is not easy to provide unique 
theoretical reasoning regarding the shapes of the halos. The situation gets 
more complicated if we include the fact that all halos will not all have the 
same mass, core radius etc and we have to modify our equations by integrating 
over the abundance of halos with a given value of mass, core radius etc. This 
brings in more ambiguities and depending on the assumptions we make for each 
of these components [e.g, abundance for halos of a particular mass could be 
based on Press-Schecter or Peaks formalism], and the final results have no real 
significance. It is, therefore, better [and probably easier] to attack the question 
based on the evolution equation for the correlation function rather than from 
"physical" arguments for density profiles. This is what we shall do next. 

4. Self-similar evolution 

Since the above discussion motivates us to look for correlation functions of 
the form ^(a, x) = a'^L{x), we will start by asking a more general question: Does 
equation (^) possess self-similar solutions of the form 




X 
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where q = xa " ?. Defining Q = \nq = X — a A and changing independent 
variables to from {A, X) to {A, Q) we can tranform our equation (j^) to the form: 

Using the relations {d^/dA)Q = (3^, {d^/dQ)A = {1/ F){dF/dQ) we can rewrite 
this equation as 



ci + Ki)]^ FdQ 
The right hand side of this equation depends only on Q and hence will vanish if 
differentiated with respect to A at constant Q. Imposing this condition on the 
left hand side and noticing that it is a function of ^(a,x) we get 

( II ) -i(Left Hand Side) = (31) 

To satisfy this condition we either need (i) {d^/dA)Q = = implying /3 = 
or (ii) the left hand side must be a constant. Let us consider the two cases 
separately. 

(i) The simpler case corresponds to /5 = which implies that ^(a, x) = F{Q). 
Setting /? = in equation (|30|) we get 



d^\_ 3(1 + 0M0 (32) 



\dQ) [a + h{i) 

which can be integrated in a straightforward manner to give a relation between 
q = exp Q and ^ : 

= go(i + 0''/'v(0-°/' 



Given the form of h[^{a,x)], this equation can be in principle inverted to 
determine ^ as a function of g = xa~". 
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To understand when such a solution will exist, we should look at the limit 
of ^ <^ 1. In this limit, when linear theory is valid, we know that h ^ (2/3)^ 
[see [Peebles 1980|| . Using this in equation (|33|) we get the solution to be 



In^ = — (2/a;) Ing or 

^ oc g-^ oc x-^a^ oc a'^x'^''-^^'^ (33) 

with the definition a = 2/{n + 3). This clearly shows that our solution is valid, 
if and only if the linear correlation function is a scale-free power law. In this 
case, of course, it is well known that solutions of the type ^(a,x) = F{q) with 
q = xa exists. [Equation (^) gives the explicit form of the function F{q)]. 

This result shows that this is the only possibility. It should be noted that, 
even though we have no explicit length scale in the problem, the function C,{q) 
— determined by the above equation — does exhibit different behaviour at 
different scales of nonlinearity. Roughly speaking, the three regimes in equation 
( pTD translates into nonlinear density contrasts in the ranges 6 < 1,1 < 6 < 200 
and 6 > 200 and the function ^(g) has different characteristics in these three 
regimes. This shows that gravity can intrinsically select out a density contrast 
of (5 ~ 200 which, of course, is well-known from the study of spherical tophat 
collapse. 

(ii) Let us next consider the second possibility, viz. that the left hand side 
of equation (RDI) is a constant. If the constant is denoted by fi, then we get 



F = Faqf" and 

P^-3{l + 0h{0=f^(y^+f^h^ (34) 



which can be rearranged to give 

{f3 - a/x)^ 



(35) 



3 + (/i + 3)e 

This relation shows that solutions of the form ^{a,x) = F{x/a'^) with (3 



- 19 - 



is possible only if h[^{a,x)] has a very specific form given by (|35|). In this form, 
h is a monotonically increasing function of ^(a, x). There is, however, firm 
theoretical and numerical evidence ( [Hamilton et al. 1991| ; Padmanabhan 1996|) 
to suggest that h increases with ^{a,x) first, reaches a maximum and then 
decreases. In other words, the h for actual gravitational clustering is not in the 
form suggested by equation (P3|). We, therefore, conclude that solutions of the 
form in equation ( |^^ with (3^0 cannot exist in gravitational clustering. 

By a similar analysis, we can prove a stronger result: There are no solutions 
of the form ^(a, x) = ^{x/F{a)) except when F{a) oc a". So self-similar 
evolution in clustering is a very special situation. 

This result, incidentally, has an important implication. It shows that 
power-law initial conditions are very special in gravitational clustering and may 
not represent generic behaviour. This is because, for power laws, we have a 
strong constraint that the correlations etc can only depend on g = xa"^/*-"'"'"^^ 
For more realistic — non-power law — initial conditions the shape can be 
distorted in a generic way during evolution. 

All the discussion so far was related to finding exact scaling solutions. It is 
however possible to find approximate scaling solutions which are of practical 
interest. Note that we normally expect constants like a, f3, /i etc to be of order 
unity while ^ (a, x) can take arbitrarily large values. If ^ (a, x) ^ 1 then equation 
( ^5l) shows that h is approximately a constant with h = {(3 — afi)/{fi + 3). In 
this case 

e(a, x) = aPF{q) oc a^q^" oc a^^-^'^^a;^ oc d'^^+'^^x^' (36) 

which has the form ^(a, x) = a^^F{a^x) which was obtained earlier by directly 
integrating equation (|^) with constant h. We shall say more about such 
approximate solutions in the next section. 
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5. Units of the nonlinear universe 

Having reached the conclusion that exact solutions of the form 
^(a, x) = a^G{x) are not possible, we will ask the question: Are there 
such approximate solutions ? And if so, how do they look like ? We will see 
that such profiles — which we shall call "pseudo-linear profiles" — that evolve 
very close to the the above form indeed exist. In order to obtain such a solution 
and check its validity, it is better to use the results of section 2.1 and proceed as 
follows: 

We are trying to find an approximate solution of the form ^(a,x) = a^G(x) 
to equation (^. Since the linear correlation function ^i(a,x) does grow as 
at fixed x, continuity demands that ^{a,x) = .^^(a, x) for all a and x. [This 
can be proved more formally as follows: Let ^ = a^G{x) and = a'^Gi{x) 
for some range Xi < x < X2- Consider a sufficiently early epoch a = at 
which all the scales in the range {xi^X2) are described by linear theory so that 
^(aj,x) = ^L(aj,x). It follows that Gi{x) = G{x) for all xi < x < X2- Hence 
^{a,x) = ^2,(a,x) for all a in xi < x < X2. By choosing a, sufficiently small, we 
can cover any range (xi,X2). So ^ = for any arbitrary range. QED]. Since we 
have a formal relation ([T0|) between nonlinear and linear correlation functions, 
we should be able to determine the form of G{x). 

To do this we shall invert the form of the linear correlation function 
^^(a,/) = a?G{l) and write / = G^^{a^'^S,L) = F{a''^S,L) where F is the inverse 
function of G. We also know that the linear correlation function ^^(a, /) at scale 
/ can be expressed as V[^(a,x)] in terms of the true correlation function ^(a, x) 
at scale x where 

/ = x(l+e(a,x))i/3 (37) 



So we can write 



/ = F 
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ra(a,/)i 


= F 


rV[e(x,a)]l 









(3^ 



But X can be expressed as x = F[^i(a, x)/a^]; Substituting this in ( p?] ) we have 



/ = F 



1 + e 



1/3 



(39) 



From our assumption C,l{0',x) = C,{ci,x) ; therefore this relation can also be 
written as 

A 1/3 



/ = F 



(40) 



Equating the expressions for / in (^Hj) and pO|) we get an implicit functional 
equation for F: 

1/3 





= F 











which can be rewritten as 



F 




F 







:i + 



1/3 



(41) 



(42) 



This equation should be satisfied by the function F if we need to maintain the 
relation ,^(a, x) = .^^(a, x). 

To see what this implies, note that the left hand side should not vary with a 
at fixed ^. This is possible only if F is a power law: 



F(o = Ar 

which in turn constrains the form of V(^) to be 

v(0 = e(i + 0'/'"^ 



(43) 



(44) 



Knowing the particular form for V we can compute the corresponding h{^) from 
the relation 

^InV 9 1 

(45) 



3 (1 + MO 
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For the V(^) considered in equation (|^) we get 

h = , , (46) 

3 + (3 + l/m)e ^ ' 

which is the same result obtained by putting /? = 2 , a = in equation (pSl). We 

thus recover our old result — as we should — that exact solutions of the form 

i{a,x) = ^l{o,,x) = a? G{x) are not possible because the correct V(^) and h{S,) 

do not have the forms in equations and (^) respectively. But, as in the 

last section, we can look for approximate solutions. 

We note from equation (^31) that for ^ ^ 1, we have 

V(0 = ^(1+1/3™); F(0 oc G{1) oc e'/" (47) 

This can be rewritten as 

V(0 = f ; F{i) oc i^l'^^'^-^)- G(0 oc f /("--i) (48) 

In other words if V(^) can be approximated as we have an approximate 
solution of the form 

e(a, x) = G{x) = x^^''-^) (49) 
Since the V in equation (|12|) is well approximated by the power laws in (0) so 



that 

V(0 oc (1 ^ ^ ^ 200) (50) 

oc f (200 < (51) 

we can take u = 1/3 in the intermediate regime and = 2/3 in the nonlinear 
regime. It follows from (HSI) that the approximate solution should have the form 



F(0 cx ^ (1 < e < 200) (52) 
oc i (200 < (53) 
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This gives the approximate form of a pseudo-hnear profile which will grow as o? 
at all scales. 

There is another way of looking at this solution which is probably more 
physical and throws light on the scalings of pseudo-linear profiles. We recall 
that, in the study of finite gravitating systems made of point particles and 
interacting via newtonian gravity, isothermal spheres play an important role. 
They can be shown to be the local maxima of entropy [ see [Padmanabhan 1990 1 
and hence dynamical evolution drives the system towards an (1/x^) profile. 
Since one expects similar considerations to hold at small scales, during the late 
stages of evolution of the universe, we may hope that isothermal spheres with 

profile may still play a role in the late stages of evolution of clustering in 
an expanding background. However, while converting the profile to correlation, 
we have to take note of the issues discussed in section 2. In the intermediate 
regime, dominated by scale invariant radial collapse ( Padmanabhan 1996| ), the 
density will scale as the correlation function and we will have ^ oc (l/s^). On 
the other hand, in the nonlinear end, we have the relation 7 = 2e — 3 [see 
equation ( pi]) ] which gives oc (1/2;) for e = 2. Thus, if isothermal spheres 
are the generic contributors, then we expect the correlation function to vary as 
(1/x) and nonlinear scales, steepening to (l/x^) at intermediate scales. Further, 
since isothermal spheres are local maxima of entropy, a configuration like this 
should remain undistorted for a long duration. This argument suggests that a ^ 
which goes as at small scales and (l/x^) at intermediate scales is likely to 

be a candidate for pseudo-linear profile. And we found that this is indeed the 
case. 

To go from the scalings in two limits given by equation (|5^) to an actual 
profile, we can use some fitting function. By making the fitting function 
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sufficiently complicated, we can make the pseudo-linear profile more exact. We 
shall, however, choose the simplest interpolation between the two limits and try 
the ansatz: 

F(z) = ^ , ^ — (54) 

where A and B are constants. Using the original definition / = F\^i^ja^\ and the 
condition that f = we get 

^ I (55) 



ila^{^ila^^B) 

This relation implicitly fixes our pseudo-linear profile. Solving for ^, we get 

(56) 

with L = AA/B'^. Since this profile is not a pure power law, this will satisfy the 
equation (|3) only approximately. We choose B such that the relation 



is satisfied to greatest accuracy at a = 1. 

This approximate profile works reasonably well. Figures 1 and 2 show this 
result. In figure 1 we have plotted the ratio F(V(0/a^)/i^(^/a^) on the x-axis 
and the function (1 + ^)^/^ on the y-axis. If the function in (^) satifies equation 
( ^2]) exactly, we should get a 45-degree line in the figure which is shown by a 
dashed line. The fact that our curve is pretty close to this line shows that the 
ansatz in (|56|) satisfies equation ( ^2|) fairly well. The optimum value of B chosen 
for this figure is -B = 38.6. When a is varied from 1 to 10^, the percentage of 
error between the 45-degree line and our curve is less than about 20 percent in 
the worst case. It is clear that our profile in (p^ satisfies equation ( pT]) quite 
well for a dynamic range of 10^ in a^. 
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2 4 



F(V(f)/a")/F(f/a2) 

Fig. 1. — The approximate solution to the functional equation determining the pseudo-linear 
profile is plotted. See text for discussion. 

Figure 2 shows this result more directly. We evolve the pseudo linear profile 
form = 1 to 1000 using the NSR, and plot [^{a^x) / a?] against x. The 
dot-dashed, dashed and two solid curves (upper one for = 100 and lower 
one for = 900) are for = 1,9, 100 and 900 respectively. The overlap of 
the curves show that the profile does grow approximately as a^. Also shown 
are lines of slope —1 (dotted) and —2 (solid); clearly ^ oc x^^ for small x and 
^ oc x"^ in the intermediate regime. 

We emphasis that we have chosen in equation (|56|) the simplest kind of 
ansatz combining the two regimes and we have used only two parameters A and 
B. It is quite possible to come up with more elaborate fitting functions which 
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Fig. 2. — The dot-dashed, dashed and two sohd curves (upper one for — 100 and lower 
one for €? — 900) are for a? — 1, 9, 100 and 900. The dotted straight hne is of slope -1 and 
the solid one is of slope-2 showing both the Xjx and 1 / regions of the profile 



-27- 



will solve our functional equation far more accurately but we have not bothered 



to do so for two reasons: (i) Firstly, the fitting functions in equation (|TT]) for 
V{z) itself is approximate and is probably accurate only at 10-20 percent level. 
There has also been repeated claims in literature that these functions have 
weaker dependence on n which we have ignored for simplicity in this paper, (ii) 
Secondly, one must remember that only those ^ which correspond to positive 
definite P{k) are physically meaningful. This happens to be the case our 
choice [which can be verified by explicit numerical integration with a cutoff at 
large x] but this may not be true for arbitrarily complicated fitting functions. 
Incidentally, another simple fitting function for the pseudo-linear profile is 

A' 

^{a,x) = (58) 
(x/L')[(x/L') + 1] 

with A! = and L' = L/A. 

If a more accurate fitting is required, one can obtain it more directly from 
equation (|16D . Setting = 2 in that equation predicts the instantaneous spatial 
slope of ^(a, x) to be 

ain^(a,x) _ 2 3(1 I _ ^ ) (59) 

dhix h[S,{a,x)] i{.o,^x) 

which can be integrated to give 

ln^= f^"'"' 

^ L Ji[L] ^(2 -3h)-3h ^ ' 

at a = 1 with L being an arbitratry integration constant. Numerical integration 
of this equation will give a profile which is varies as (1/x) at small scales and goes 
over to (l/x^) and then to (l/x^), (1/x^).... etc with an asymptotic logarithmic 
dependence. In the regime ^(a,x) > 1, this will give results reasonably close to 
our fitting function. 
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It should be noted that equation (^21) reduces to an identity for any F, in 
the hmit ^ — since, in this hmit V{z) ^ z. This shows that we are free to 
modify our pseudo-hnear profile at large scales into any other form [essentially 
determined by the input linear power spectrum] without affecting any of our 
conclusions. 

Finally, we will discuss a different way of thinking about pseudolinear profiles 
which may be useful. 

In studying the evolution of the density contrast (5(a,x), it is conventional to 
expand in in term of the plane wave modes as 

6{a, x) = ^(g, k) exp(zk ■ x) (61) 

k 

In that case, the exact equation governing the evolution of 6 {a, k) is given by 
( Peebles 198"0D 



da? 2a da 2a^ 
where A denotes the terms responsible for the nonlinear coupling between 
different modes. The expansion in equation (|6TD is, of course, motivated by the 
fact that in the linear regime we can ignore A and each of the modes evolve 
independently. For the same reason, this expansion is not of much value in the 
highly nonlinear regime. 

This prompts one to ask the question: Is it possible to choose some other set 
of basis functions Q{a, x), instead of exp ik ■ x, and expand 6{a, x) in the form 

6{a,x.) = ^6a{a) Q{a,x.) (63) 

a 

SO that the nonlinear effects are minimised ? Here a stands for a set of 
parameters describing the basis functions. This question is extremely difficult 
to answer, partly because it is ill-posed. To make any progress, we have to 
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first give meaning to the concept of "minimising the effects of nonhnearity" . 
One possible approach we would like to suggest is the following: We know 
that when 5(a, x) <^ l,then 5(a,x) oc a F{x.) for any arbitrary -F(x); that is 
all power spectra grow as in the linear regime. In the intermediate and 
nonlinear regimes, no such general statement can be made. But it is conceivable 
that there exists certain special power spectra for which P(k, a) grows (at 
least approximately) as even in the nonlinear regime. For such a spectrum, 
the left hand side of ( 162|) vanishes (approximately); hence the right hand side 
should also vanish. Clearly, such power spectra are affected least by nonlinear 
effects. Instead of looking for such a special P(/c, a) we can, equivalently look 
for a particular form of ^(a,a;) which evolves as closely to the linear theory as 
possible. Such correlation functions and corresponding power spectra [which are 
the pseudo-linear profiles] must be capable of capturing most of the essence of 
nonlinear dynamics. In this sense, we can think of our pseudo-linear profiles as 
the basic building blocks of the nonlinear universe. The fact that the correlation 
function is closely related to isothermal spheres, indicates a connection between 
local gravitational dynamics and large scale gravitational clustering. 

6. Conclusions 

It seems reasonable to hope that the late stage evolution of collisionless 
point particles, interacting via newtonian gravity in an expanding background, 
should be understandable in terms of a simple paradigm. This paper [ as the 
title implies! ] tries to realise this dream within some well defined framework. 
It should be viewed as a tentative first step in a new direction which seems 
promising. 

There are three key points which emerge from this analysis. The first is 
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the fact that we have been able to find approximate correlation functions 
which evolve preserving their shapes. We achieved this by looking at the 
structure of an exact equation which obeys certain nonlinear scaling relations. 
As we emphasised before, the existence of such special class of solutions to the 
equations of gravitational dynamics is an important feature. 

Secondly, we should take note of the role played by the "isothermal" profile 
(1/x^) in our solution. Such a profile can lead to correlation functions which go 
as (1/x) at small scales and in the intermediate scales. If this profile is 

indeed "special" then one expects it to lead to a pseudo-hnear profile for the 
correlation function. Our analysis shows that there is indeed good evidence 
for this feature. If one accepts this evidence, then the next level of enquiry 
would be to ask why profiles are "special". In the statistical mechanics 

of gravitating systems, one can show that these profiles arise as end stages of 
violent relaxation which operates at dynamical time scales. Whether a similar 
reasoning holds in an expanding background, independent of the index for 
power spectrum, is open to question. This is an important issue and we hope to 
address it fully in a future work. We emphasise that our equations, along with 
NSR, naturally lead to a pseudo-linear profile, which can be interpreted and 
understood in terms of isothermal density profiles for halos; we did not have to 
assume anything a priori regarding the halo profiles. 

In a more pragmatic way, one can understand the pseudo-linear profile from 
the dependence of the rate of growth of the correlation function on the local 
slope. The NSR suggest that ^ grows (approximately) as a^/^^^J^^^^^ in the 
intermediate regime and as a^/^^<'ff+^^ in the nonlinear regime. This scaling 
shows that Ugff — —1 grows as in the intermediate regime and Ueff — —2 
grows as in the nonlinear regime. This is precisely the form our pseudo-linear 
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profile has. Also, in the intermediate regime, the correlation grows faster than 
a? if Ueff < —1 and slower than a? if n^ff > —1. The net effect is, of course, to 
straighten out a curved correlation and drive it to = — 1. Similar effect drives 
the correlations to n = —2 in the nonlinear regime, [see Pagla fc Padmanabhan 



1997 for a more detailed discussion of this aspect in the intermediate regime] 



Of course, one still needs to understand the dependence of growth rate on the 
Ueff from more physical considerations to get the complete picture. We have not 
addressed in this paper, what is the timescale over which clustering can lead to 
the psuedo-linear profile even granting that it does. This requires further study. 

The last aspect has to do with what one can achieve using the pseudo-linear 
profiles. In principle, one would like to build the nonlinear density field through 
a superposition of pseudo-linear profiles but this is a mathematically complex 
problem. As a first step one should understand why the nonlinear term in 
equation (^) is subdominant for such a profile. This itself is complicated since 
we have only fixed the power spectrum — but not the phases of the density 
modes — while the nonlinear terms do depend on the phase. Again, we hope to 
investigate this issue further in a future work. 
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